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We conduct a numerical study of the dynamical behavior of a system of three-dimensional 
"crosses", particles that consist of three mutually perpendicular line segments of length cr rigidly 
joined at their midpoints. In an earlier study [W. van Ketel et at, Phys. Rev. Lett. 94, 135703 
(2005)] we showed that this model has the structural properties of an ideal gas, yet the dynamical 
properties of a strong glass former. In the present paper we report an extensive study of the dynam- 
ical heterogeneities that appear in this system in the regime where glassy behavior sets in. On the 
one hand, we find that the propensity of a particle to diffuse is determined by the structure of its 
local environment. The local density around mobile particles is significantly less than the average 
density, but there is little clustering of mobile particles, and the clusters observed tend to be small. 
On the other hand, dynamical susceptibility results indicate that a large dynamical length scale 
develops even at moderate densities. This suggests that propensity and other mobility measures are 
an incomplete measure of dynamical length scales in this system. 

PACS numbers: 61.43.Fs, 61.20.Lc, 05.20.Jj 



I. INTRODUCTION 

There exist a bewildering variety of theories for the 
ffkss transition (see e.g. [Ill i, i, i, H, 0, i, i, M, HH, 
Roughly speaking, one can distinguish between 
two main classes. Theories belonging to the first class 
are based on the assumption that static, structural cor- 
relations in the fluid are ultimately responsible for the 
occurrence of structural arrest. Theories that belong to 
the second class assume that purely kinetic factors con- 
trol the onset of glassy behavior. It is probably fruitless 
to search for the "true" theory of the glass transition, 
because not all experimental glasses appear to be equiv- 
alent [13, [l^. However, it is important to disentangle, 
as much as possible, the roles of structural correlations 
and of purely kinetic effects in the absence of such corre- 
lations. 

Recently, we reported simulations that provided evi- 
dence that it is possible to observe glassy behavior in 
a model s yste m that has the structural properties of an 
ideal gas [ia|. As the particles in an ideal gas have no 
static structural correlations, dynamical arrest in this 
system is a purely kinetic effect. The model system 
we explore consists of particles made of three mutually 
perpendicular line segments of length cr, rigidly joined 
at their midpoints. These three-dimensional "crosses" 
generalize the hard-needle model developed to study 
topological effects on rotational and translational diffu- 
sion [rn . [ist , as has already been implicitly [1] or explic- 
itly [l0| suggested. A lattice-based version of the hard- 
needle system has already been studied by several groups 
as a model for orientational glass formers [l^, dl], [22, [l^] ■ 
Renner et al. [l^l simulated line segments that can ro- 
tate around fixed lattice points. The system enters a 
non-ergodic glassy phase at finite segment length, but 
since it has an ideal static behavior the standard mode- 



coupling theory (MCT) of the glass transition is inap- 
plicable. However, an extension of MCT that includes 
torque-torque contributions doespredict a glass transi- 
tion for these lattice rotators [13, [l^. Closer to the 
present model is the thin line segments with fixed but 
random orientations, whose dynamics was studied by 
Szamel et al. (2^ . [25| . Using a mean- field approximation, 
they found that the transverse motion of the line seg- 
ments decreases severely with increasing segment length, 
because of entanglement (tube constraints), while the 
motion along the orientation of the lines is not affected 
by such constraint. 

Since the crosses have zero volume and thus zero ex- 
cluded volume, all static thermodynamic quantities are 
exactly known. By random insertion one can trivially 
generate a representative equilibrium configuration at 
any density. As our model is an ideal gas, the onset 
of glassy behavior takes place within a single thermody- 
namically stable phase. Hence we need not worry that 
the dynamics in the glassy phase be obscured by the slow 
nucleation of another phase. We can also safely ignore 
the "Kauzmann paradox" [2^ , which states that the glass 
transition takes place when the entropy of the fluid phase 
threatens to drop below that of the crystal phase. The 
present model has no crystal phase nor for that matter 
any ordered phase. Nonetheless, its dynamics is highly 
nontrivial. 

The rest of the paper is organized as follows. In Sec.HIl 
we describe the model and the simulation algorithm. 
Collision as well as diffusion properties are presented in 
Sec. nil Al and the self-intermediate scattering function in 
Sec. HUB] We investigate in details the effect of the local 
environment on the mobility of particles and clustering of 
the "mobile" particles in Sec. IIII Cl while Sec. |TTTD] con- 
cerns itself with the density and wave- vector dependence 
of the four-point susceptibility. Finally, we conclude in 
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Sec. IIVI with a summary of the important findings. 



II. SIMULATION TECHNIQUE 

We simulate a three-dimensional system at constant 
number of particles N, volume V and temperature T un- 
der Newtonian dynamics. The particles consist of three 
mutually perpendicular line segments of length a rigidly 
joined at their midpoints. We choose the initial center 
of mass positions of the particles at random in the cubic 
simulation box. The box volume V — Na^/p is set by 
the choice of N and the number density p. The cross ori- 
entations are also randomly distributed. For numerical 
convenience we reject configurations having two crosses 
with almost identical orientations to within an angle of 
lO"** radians. Assuming truly random orientations, the 
probability of having such closely aligned pairs of crosses 
is less than one part in 10^ for the system sizes consid- 
ered. Hence, the effect of this choice should be negli- 
gible. We also neglect rotational motion, which corre- 
spond to having crosses with an infinite moment of iner- 
tia, so they preserve their initial orientation throughout 
the simulation. This allows us to analytically compute 
the time before the next collision, which leads to large 
computational efficiency gains. The initial velocities are 
randomly drawn from a Maxwell-Boltzmann distribution 
and shifted to set the center of mass velocity to zero. We 
choose a as the unit of length, the thermal energy ksT 
as the unit of energy, and the particle mass m as the 
unit of mass. This results in time t being expressed in 
units of {ksT /ma'^)~^^'^ . Simple periodic boundary con- 
ditions are used in all three directions. The dynamical 
rules are simple: between collisions, the particles move 
ballistically, while when two line segments collide, the 
component of relative velocity perpendicular to the plane 
of the two line segments is reversed. 

We use an event-driven algorithm [2^ , wherein future 
collision events are stored in a binary tree structure and 
the particle positions arc updated asynchronously to the 
time of the next collision event. Because of the extreme 
anisotropy of the crosses, both spherical neighbor lists 
and cell structures are inefficient at high densities. In- 
stead we consider spherocylinders around each line seg- 
ment and create neighbor lists from the spherocylinder 
overlaps. "Events" in our algorithm are not only colli- 
sions, but also neighbor list updates. These occur when- 
ever the center of mass of a particular cross has moved 
by more than half the spherocylinder radius since that 
cross's neighbor list was last created. To limit the search 
for spherocylinder overlaps while creating the neighbor 
lists, we consider a cubic cell structure based on the 
center of mass of the crosses. We limit the size of the 
event tree by setting a time (typically five times what it 
would take a particle with the average speed to ballisti- 
cally cross the neighbor list cutoff length) beyond which 
events are not entered in the tree structure, which also 
sets the longest survival time of a neighbor list. This 



ensures that if a particular cross does not undergo any 
collision within this interval we still correctly identify fu- 
ture events that involves it. 

For the very rare case of quasi-simultaneous collisions, 
the behavior of the program is unpredictable. Depend- 
ing on the exact sequence of instructions, a future event 
can behave like a past event and vice versa. We avoid 
this problem by discarding events that are separated by 
less than 10^^^ time units from a previous event. Since 
this time is much smaller than the average time between 
collisions even at the highest density considered in this 
study, this artificial exclusion does not affect the statis- 
tical analysis of our data. 

After a collision, all events involving the colliding pairs 
are removed from the event tree and new future events 
are generated from their respective neighbor lists. When 
the event is a neighbor list update, all events involving 
this particle are removed and the list is recreated anew. 
When the next event is later than the time at which we 
are required to calculate any property of the system, we 
update the positions and velocities of all the particles 
to that time without changing the event list, since by 
definition the next event is later than this time. 

For the highest densities and largest system sizes con- 
sidered in this work, the number of collision events in a 
single run often exceeds 10^°. On a 2 GHz AMD Opteron 
linux desktop using an Intel Fortran compiler, the CPU 
time required for 10^° collisions to take place in a sys- 
tem of 4096 crosses at p = 20 is about 25 hours. The 
two most costly operations are finding the future colli- 
sions and filling the spherocylindrical neighbor list. The 
optimum performance is observed when the (density de- 
pendent) spherocylinder radius is chosen such that the 
average number of neighbors is about 30. For a smaller 
radius the neighbor list is updated more often, while for 
a larger radius future collisions are found among a larger 
set of possible interactions. 

A random insertion procedure gives an equilibrated 
configuration for the ideal gas, since the radial distri- 
bution function g{r) is flat. But the dynamics retains a 
long memory and the structural relaxation slows down 
exponentially with density. For this reason it is more 
efficient to perform the averaging by choosing statically 
independent starting configurations and running them on 
different cores. For most of our simulations we use 512, 
1728, and 4096 particles with p varying from 1 to 30. 
All the simulations up to p = 20 and N — 4096 are 
run for at least 10^ collision times or until the small- 
est nonzero wave vector q = 2'kV~^/'^ component of the 
dynamic structure factor S{q,t) / S{q,Q) has decayed to 
1/e, whichever is smaller. The runs with p > 20 and 
N > 4096 are not sufficiently long to satisfy the second 
condition, so they are only used to determine quantities 
measured on shorter time or length scales. The averag- 
ing procedure employed still guarantees the validity of 
these results. Finite-size effects are found to be negligi- 
ble for all static and two-point quantities in the density 
regime under study, but four-point correlations exhibit 
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FIG. 1: (Color online) Average collision time Tcoi from simu- 
lation (points) compared to the the kinetic theory prediction 
(dahsed line) . Inset : Collapse of the velocity autocorrelation 
function Z{t) after rescaling time by Tcoi- 



sizeable size dependence. This will be discussed further 
in Sec. HITD] 




FIG. 2: (Color online) Time evolution of the MSD for p= 1, 2, 
5, 7, 10, 12, 15, 17, and 20, from left to right. Superimposed to 
the long time part of p = 20 is a linear fit whose slope is used 
to calculate the diffusion coefficient. The error is smaller than 
the symbol size. Inset : The diffusion coefficient decreases 
with density exponentially. The dashed line is a fit D ~ 
exp[-Ay*p] with AV* = 0.42 for p > 5. 



III. RESULTS 



A. Collisions and Diffusion 

By construction the static properties of the system are 
those of an ideal gas. Moreover, for the density range 
considered the short-time dynamics agrees within errors 
with a mean- field kinetic theory. Fig. [T] shows in fact that 
the average time between rod collisions is indistinguish- 
able from the analytical prediction Tcoi — 4/(9p7r^/^) 
(see Appendix). After only a few collisions the par- 
ticle velocities become nearly completely uncorrelated, 
as gathered from the velocity autocorrelation function 
Z{t) = {vj{t) • Vj(0))/(|i;p) (Fig. [U inset). The nearly 
perfect collapse of Z(t) after rescaling time by t^oi m 
Fig. [1] shows this process to be rather general. The small 
negative dip of Z{t) that follows at high density is the 
caging signature and corresponds to the bouncing back of 
a particle after colliding with a neighbor. As far as struc- 
ture and short-time dynamics are concerned the system 
thus behaves rather ideally. 

On longer timescales the physics is quite different. 
Fig. [2] shows that the mean-square displacement (MSD) 
[{Ar^{t)) EE {\rj{t) - rj(O)p)] between the initial bal- 
listic regime [(Ar^(i)) ~ i^] and the diffusive regime 
[(Ar^(t)) ^ 6Dt], where D is the diffusion coefficient, 
develops a plateau for increasing densities as in super- 
cooled fluids. But contrary to structural liquids there is 
no upper limit to packing, so the transition away from 



the ballistic regime takes place at ever shrinking length 
and time scales with increasing density. Instead of con- 
verging at a single length scale set by the repulsive core, 
as is the case in structural glass formers, the crossover 
plateau thus keep lowering with the slowdown. With the 
end of the plateau region, the system enters the diffu- 
sive regime on a timescale that grows exponentially with 
density. This suggests that the rate-limiting step for dif- 
fusion is the creation of "free volume" around a particle, 
such that the topological constraints inhibiting its mo- 
tion are relieved. For an ideal gas the probability to 
open up a volume AV* by a spontaneous fluctuation is 
^ exp{—pAV*). The exponential density dependence of 
D thus suggests that a cavity with volume AV* ~ QA2a^ 
is needed to enable diffusion. This behavior is very differ- 
ent from the algebraic density dependence observed for 
the rotational diffusion in systems of tethered, rotating 
needles 20] . It is also unlike that of structural athermal 
systems such as hard spheres, where a power law is ob- 
served at modest undercooling [2^ . Exponential slowing 
down is more akin to what is obtained in strong glass 
formers. 



B. Self-intermediate scattering function 

The decay of density fluctuation on different length 
scales is best studied by the incoherent self-intermediate 
scattering function Fs{q,t) = X]j 6xp{«q • [rj(t) — 
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FIG. 3: (Color online) Fs{q,t) decay at (a) the microscopic wave vector gcagc as well as its Tq collapse at (b) q — 0.8n and (c) 
q = 2-K for p = 5, 7, 10, 12, 15, 17, and 20. The solid line is a stretched exponential fit to Aexp[—(t/Ta)'^] for Fs{q,t) between 
0.025 and 0.975 with (b) ^4=0.975, /3=0.917 and (c) y4=0.963, /3=0.662. Insets: Short time decay of Fs{q,t) with additional 
p—22, 25, and 30. High-density (p > 20) estimates for are obtained by forcing the early a decay onto the master curve from 
the low-density data. The solid line is the free-particle decay form with 2.4rcoi, as described in the text. 
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FIG. 4: (Color online) (a) Structural relaxation time extracted from Fs{q,t) for various wave vectors. The dashed line is 
an exponential fit to g = tt with exponent 0.43. (b) Exponents extracted from Fs{q,t) and X.lilt't) at p = 20, as described 
in the text of Sec. IIII Al and Sec. IIIIDl (c) Rescaling of Ta{q) by the difi'usive limit Dq^ to evaluate the transport coefficient 
decoupling. The solid line emphasizes the Fickian limit Dq^Ta{q) = 1, while the dashed lines show the small wavevector limit 
\f2Dq for the two lowest densities. 



rj(0)]}), where g = |q| as reported in Fig. [31 In standard 
glass formers this correlation function bears the signature 
of two different dynamical regimes in the microscopic re- 
laxation. On times of the order of Tcoi, ballistic motion 
gives way to the (5 plateau associated with caging; on 
longer times scales, a structural rearrangements allow a 
particle to escape the cage formed by its neighbors. The 
typical timescale (g) over which this last process takes 
place is defined as the time when Fs{q,t) has decayed 
to 1/e. Here, Tq, increases exponentially with density 
(Fig. 11^) . This supports the assumption that an infinite 
cross density is necessary to obtain complete dynamical 
arrest. The length scale at which caging and structural 
relaxation are best separated is the caging diameter. In 
structural glass formers it also corresponds to the first 
peak of the structure factor, but since the crosses do not 
exhibit any static structure,we approximate it instead by 
the average spacing between particles gcago = 27rp^/'^. 
The growing separation between the two timescales with 



density can be observed in Fig. [3^. However, in spite of 
there being a difference of over three orders of magni- 
tude between Tcoi and Tq at p = 20, the apparition of a 
plateau is still incomplete. Higher densities are necessary 
to observe a better delineated structure. 

An analysis of the structural relaxation process shows 
that when time is rescaled by Ta{q) for a fixed q, Fs{q, t) 
collapse onto a single master curve with increasing accu- 
racy as the system gets denser (Fig. |3}d-c). This time- 
density scaling of the long-time decay is highly non- 
trivial. It has been argued that the collapse of the a- 
relaxation curves is one of the outstanding characteris- 
tics of the structural glass transition that is reproduced 
by MCT 0, Q . Standard MCT being here inappUcable 
for lack of static correlations the phenomenon is clearly 
more generic. Long-time relaxation in glassy systems are 
often described by a stretched exponential Kohlrausch- 
Williams- Watts (KWW) form Fs{q,t) ~ e-^*/^")" , where 
the stretching exponent (3 is not to be confounded with 
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the /3-relaxation regime. The long length scale limit is 
properly captured at low wave vectors, as seen in Fig.|3|D. 
The stretching exponent (3 then approaches unity. This 
corresponds to an exponential decay of Fs (q, t) and is 
consistent with the diffusive dynamics of simple flu- 
ids. At microscopic length scales the KWW fit is also 
rather successful, as shown in Fig. [3)3-c, though no sin- 
gle parametrization of the functional form is suitable for 
the entire decay range [l^. In particular, the fitting 
form does not capture the long-time tail, which falls off 
faster than expected from a fit to the body of the decay. 
For the bulk of the decay however, a singular behavior 
is observed. In other model glass formers such as sil- 
ica [2§| and binary Lennard- Jones (LJ) [sO], P > 0.75 for 
wavevectors around gcagc- For the crosses the decay is 
further stretched with (3 w 0.5. This suggests that the 
structural relaxation arises from a broader characteristic- 
time distribution of relaxation processes. 

MCT further predicts that the end of the plateau 
bends down following a von Schweidler form [j, Q 

Fs{q,t) ^ Uq) ^ Bh{q){t/T^)\ (1) 

where fc{q) is the plateau height and both B and h{q) are 
independent of time. Equation [1] approaches a stretched 
exponential form in the large-q limit [3l|. For densities 
considered here this form is not obviously appropriate, 
since no convincing plateau has yet developed. This 
leaves fdq) as a free fitting parameter to extract the ex- 
ponent b from the decay shoulder at p = 20, as reported 
in Fig. UJ). Though coarse, this treatment will be useful 
when we return to this issue in Sec. IIII Dl 

A feature not part of the canonical glass analysis is the 
short-time collapse of Fs{q,t), as presented in Ref. [T6| 
and depicted in the insets of Fig. [3Jd-c. At short times 
the particles' ballistic movement leads to an initial Gaus- 
sian decay of Fg {q, t) . This regime ends when the "free" 
crosses collide with the "cage" formed by their neighbors 
at time Tcoi on average. Using 

^^•=(9, Tcoi) = exp {-kBTq\lJ2m) (2) 

and the scaling of with density, one can parameter- 
ically plot where the change of regime from ballistic to 
collisional should take place for various densities. Since 
this does not correspond directly to a particular feature 
of Fs{q,t), let's consider a larger value than Tcoi to de- 
scribe the observed change in regime. Mobile particles 
have more free space around them (see Sec. IIII Cp and 
contribute longer to the free decay of Fs {q, t) , so this is 
not unreasonable. Equation [2] with 2.4tco1 indeed cap- 
tures the regime change at early times, as seen the insets 
of Fig.[3l3-c. This time parameter is close to the first zero 
of Z{t), another metric for the onset of caging. This ex- 
planation is rather system specific, so this collapse is not 
expected to be observed in other glass-forming systems. 



C. Dynamical heterogeneity 

Various transport properties correspond to different 
moments of the distribution of microscopic times, so 
their decoupling at a particular wavevector is associated 
with the growth of dynamical heterogeneity on the cor- 
responding length scale [H, [11] . We first probe this ef- 
fect using the wavevector dependence of Ta{q) rescaled 
by Dq^ and then looking for the onset of decoupling 
Dq^Ta{q) > 1. As seen in Fig. [J];, for small wavevectors 
the Fickian limit Dq^Ta{q) = 1 is recovered, while at very 
high wavevectors the Gaussian decay of Fs{q, t) leads to a 
trivial ^/2Dq growth. The transition from one regime to 
the other takes place over microscopic sizes q > (/cage- In 
denser systems decoupling is more pronounced and takes 
place at increasing length scales. For the highest densi- 
ties the onset of decoupling suggests that particles have a 
coherent dynamics over distances as large as 4 — 5ct. This 
is similar to what is observed in binary LJ under simi- 
larly sluggish relaxation (32i . i34j , but here the number of 
particles involved is here an order of magnitude larger. 
We will come back to this issue in Sec. IIIIDl but note 
for now that since this size scale corresponds to the box 
dimension at these densities, it sets a computational up- 
per bound to the range of densities reasonably accessible 
through simulations. 

A number of simulation [H, IH, [13, [H, [11| and experi- 
mental [10, [4l[ studies of glass- forming systems also show 
a close relationship between the non-Gaussian behavior 
of particle displacements and dynamical heterogeneity. 
At high densities the dynamics of the crosses is indeed 
heterogeneous: only a small fraction of all particles is re- 
sponsible for a significant fraction of the total MSD be- 
tween the ballistic and the diffusive regimes, where the 
MSD plateaus. The probability distribution of particle 
displacements in Fig. [5^ shows a tail at high displace- 
ments for intermediate times, while at short and long 
times the distribution tends towards a Gaussian. Devia- 
tions can be quantified using higher-order cumulants, the 

simplest of which is the fourth-order a2{t) = ^^-2 ^^y^^ — 1. 
It vanishes when a distribution is truly Gaussian, but for 
the crosses and for structural glass formers it peaks more 
prominently and at longer times with increasing density. 
At p = 20 and time Tq.^ , when the non-Gaussian param- 
eter a2{t) reaches its maximum value, only 5% of the 
particles are responsible for nearly 30% of the MSD [3| ■ 

To look further in the microscopic features of this phe- 
nomenon, mobile and slow particles have to be identified. 
The distinction between the two types is not a sharp 
one and depends on the time interval under considera- 
tion. For both short and long time intervals all particles 
have a similar MSD and the labels lose their meaning 
altogether. Kob et al. define a critical value of the dis- 
placement at a given time beyond which the self part of 
the van Hove function deviates significantly from the cor- 
responding Gaussian approximation [35| . Particles that 
have a displacement larger than this critical value are 
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FIG. 5: (Color online) For a system at p — 20. (a) Displacement probability distribution with superimposed Gaussian fits 

P{r{t)) = 47rr-^(t) (27r{r^(t))/3) "^''^ exp(-3r^(t)/2{r^(t))) for t = 5A x 10~^ 0.059, 0.39, 15, 580, and 3.6 x 10^ from left to 
right. Arrows point to the excess probability for particles with large displacements. Inset: One-dimensional component of the 
displacement probability at t = 32 ~ Ta2 with a Gaussian and an exponential fit small and large amplitudes respectively, (b) 
Propensity probability distribution P((Ar^)|^^) for the same first five times. Inset: distribution of displacements for the 0.07% 
particles with the largest (open symbols) and smallest (closed symbols) propensities at t = 33. (c) Propensities at t = 33 shown 
as spheres centered around the initial particle positions. The spheres have a radius 2.5 times the magnitude of the individual 
particle propensities. The box has a side 3cr. 
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FIG. 6: (Color online) Mobility analysis for a system of 20,000 crosses at p = 20. (a) The radial distribution function g{r) is 
compared with the conditional distributions Qmair) and gmm{r). (b) Correlation of the mobile particle displacement directions, 
as described in the text, (c) Displacements of clustered mobile particles over Tq-j . The cone's base is at a mobile particle's initial 
position and the cone's height is twice its squared displacement over Ta2- Different shadings code for independent clusters. The 
box has side 10a. 



termed mobile. This distinction between the two regimes 
can be observed near Tq.^ in the inset of Fig. [5^, where the 
short-range Gaussian and long-range exponential separa- 
tion suggested in Ref. [12] captures the data reasonably 
well. From a different point of view, Shell et al. showed 
that the joint probability distribution of initial velocity 
component and displacement along the same direction 
can be fitted by the sum of two Gaussian functions at 
intermediate times [43]. These authors suggest that the 
relative weights of the Gaussian components can be used 
to estimate the fraction of particles that are respectively 
mobile and immobile on that timescale. In this study we 
find that different measures of heterogeneity yield essen- 
tially the same results near ■ For this reason we use a 
simpler prescription: the 5% of particles with maximum 
displacement at are termed mobile [38 1. 

In order to understand the physical origin of dynam- 
ical heterogeneities, it is important to gain insight in 



the factors that make a particular particle mobile. One 
possibility is that the distance over which a particle 
moves is sensitive to the initial velocity of that parti- 
cle, but Z{t) decays so rapidly that this could hardly 
be the whole story. Alternatively, the future mobility 
of a particle can be related to the detailed geometry 
of its initial local environment [3 SI, [H, 113, El]. To 
distinguish between the two we consider an "ensemble" 
of trajectories that initiate from the same starting con- 
figuration, but with different initial velocities. Simula- 
tions of such an "iso-configurational ensemble" (IC) al- 
low to determine whether the propensity for high mo- 
bility is related to the initial velocity or to the initial 
structure. If the former holds different particles are mo- 
bile from one trajectory to another, while if the latter 
holds the identity of mobile particles is correlated over 
different trajectories [3, HE S, SI] . For this we define 
the particle's propensity to diffuse at time t as the IC 
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average of the square of its displacement {Arf{t))ic = 
(|ri(i) — ri(0)p)i(7- At both short and long times the dis- 
tribution is expected to tend towards a delta function, 
because all starting positions are equivalent. Before a col- 
lision takes place all heterogeneities are kinetic, while for 
t 3> all possible environments are sampled. If there 
is a structural contribution to dynamical heterogeneity it 
should thus appear at intermediate times. We run repli- 
cates of identical starting configurations at p = 20 to 
look at the distribution of propensities. Figure [SJa shows 
that removing the spread due to kinetic effects indeed 
gives a thinner propensity distribution than the full dis- 
placement distribution of Fig. [5^. But the relative width 
of the displacement distribution still grows until t ^ Tq^ 
and decreases afterwards. There is thus a structural com- 
ponent to dynamical heterogeneity in the cross system. 
However, though some particles have a propensity much 
higher than others, no feature of the distribution allows 
for a separation between propensity regimes, contrary to 
structural glass formers [4J| . To see if there is nonetheless 
speciation, we look at the displacement distribution for 
the extremes of propensity. Yet in spite of having an av- 
erage propensity an order of magnitude apart, their dis- 
placement distributions at i ~ still overlap (Fig. ]Ejp 
inset). Thus only a probabilistic propensity categoriza- 
tion is possible at the particle level. But as for structural 
glass formers, it could still indicate that certain regions 
of space are structurally more mobile than others [49| . 
We consider this option in Fig. [5fc, where the spatial dis- 
tribution of particle propensities at is depicted as 
spheres centered around the initial particle position. It 
is hard to properly assess the regions of higher mobility 
directly from this representation. Though there appears 
to be some mobile "domains" , where the most highly mo- 
bile particles can be found, these are not very large and 
for the rest the mobile particles appear more or less uni- 
formly distributed over the system. This is significantly 
different from the large regions of similar propen sity that 
are observed in structural glass formers |4^. |45|. Hy, \^ . 
Either dynamically heterogeneous regions are here much 
smaller or propensity is an insufficient microscopic ob- 
servable to capture their essence in crosses. 

We nonetheless examine quantitatively possible spatial 
correlations among mobile particles with eight instances 
of a system of 20, 000 crosses dX p — 20. The radial distri- 
bution function distinguishing the mobile particles from 
the rest is compared to the featureless system- wide g{r) 
in Fig. [5^. The conditional probability of finding any 
particle at a distance r given that a mobile particle is lo- 
cated at the origin gmaif) shows a depression near r = 0. 
This indicates that mobile particles tend to be found in 
local low-density regions, as suggested by the relaxation 
mechanism presented in Sec. IIII Al The radial distribu- 
tion function of mobile particles alone gmm{f) indicates 
that they are also spatially correlated. It appears from 
this that mobile particles do organize in clusters over ex- 
tended volumes. A different measure of correlations in 
the mobile particle distribution considers the displace- 



ment directions of mobile particles. For this we define a 
correlation function 



Om{r) 



(Ar^(O) ■ Ar^(r)) 
(|Ar„J2) 



(3) 



where Arm(O) is the displacement over the time interval 
of a mobile particle considered to be at the origin and 
Ar^ (r) is the displacement of mobile particles in a spher- 
ical shell of radius r. Without correlations among the 
displacement direction of mobile particles Omir) would 
be zero, while a non-zero value indicates some degree of 
assistance between mobile particles. Fig.Hb shows a pos- 
itive Orn at small r, so mobile particles' movements are 
only correlated when they are sufficiently close together 
to be "entangled". The negative dip that follows might 
be due to poor statistics, but this cannot be resolved 
here. 

In structural glass formers mobile particles are some- 
times found in clusters with a ramified morphology ^3d\. 
The analysis done so far leaves open the possibility of 
chain-like movements for the cross model, which incites 
us to look directly at the spatial distribution of mobile 
particle clusters. Here, two mobile particles belong to a 
same "cluster" if their separation is less than a/2 in all 
directions at both initial and final times. This threshold 
is similar to the decay length scale of 5m„i(r). It ensures 
that members of a cluster share collision history over the 
entire time interval during which displacement is consid- 
ered. Most mobile particles do not belong to such a clus- 
ter and only 10% of them belong to clusters of size six or 
more; the largest cluster identified contains 14 particles. 
Figure [Ht shows clusters of six or more mobile particles 
as cones with a base centered around the particles' initial 
position and oriented along their displacement. We find 
no indication of non-compact or linear chains of mobile 
particles contrary to what was observed in simulations of 
the binary L J glass former [s^ [s^ . This allows to con- 
clude that high- mobility clusters do indeed exist and that 
they arc not only small, but also compact. But at such 
high density, though the system has undergone a signif- 
icant dynamical slowdown, collectively relaxing regions 
remain of limited spatial extent. 



D. Dynamical susceptibility 

A particularly useful quantity to discriminate between 
different models of dynamical arrest and to provide fur- 
ther information about the relaxation mechanism is the 
four-point density correlator ^ [U, M, [H, HI HI, HI 

G4(r,i) = (Ap(0,0)Ap(0,t)Ap(r,0)Ap(r,t)) 

-(Ap(0, 0)Ap(0, t)) (Ap(r, 0)Ap(r, i)), (4) 

where Ap(r, t) denotes a density fiuctuation at position r 
and time t. Gi probes the spatial correlation in the decay 
of density fluctuations at different times. The volume 
integral of G'4(r,t) is its associated susceptibility Xi{t)i 
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FIG. 7: (Color online) (a) Determination of Xii^^'t) with (TV) = 8192 at p = 5 by two diflterent approaches: direct simulation 
(full symbols) and through Eq. [8] (empty symbols), (b) System-size dependence of xT at p = 5. Inset: X4 for {^) ~ 8192 at 
constant p (squares) and constant /x (circles), (c) Dynamic susceptibility X4(9i*) for P = 20. Inset: The solid line follows 
the ballistic behavior at t < Tcoi and the dashed one the intermediate power-law regime. 



which is also a measure of the variance of the correlation 
function (Ap(0, 0)Ap(0, t)). Numerical simulations show 
that the information contained in this reduced dynamic 
susceptibility is very similar to the full four-point density 
correlator |52l |. In practice it is convenient to compute a 
phase-space correlator in terms of the self-intermediate 
scattering function 



1 ^ 



ri(t)-rj(0)] 



(5) 



From this definition we recognize that Fa{q,t) = 
(/^(q, i)). In athermal systems, the corresponding dy- 
namic susceptibility is then 



xl{q.t)=N (/,(q,i)2) -(/,(q,t))^ 



(6) 



at constant density. We use the p label, because unlike for 
the two or the full four-point correlators the susce ptib ility 

tends on the choice of simulation ensemble [541 . ISSl . 
The "true" susceptibility is obtained by keeping the 
chemical potential p, fixed instead. This can be done 
directly or using the derivative of the two-point function 



dF,{q,t) 
d\np 



(7) 



where kt is the isothermal compressibility and p refers to 
the constant chemical potential. For crosses with ksT = 
1 it reduces to 



dFsiq,t) 
dlnp 



(8) 



This result is tested for (N) = 8192 at p = 5 in Fig.lTK. 
We use the scheme described in the Appendix of Ref. [53 
on the one hand and numerical differentiation of the two- 
point function added to Xiily't) on the other. The two 
approaches agree with each other within numerical uncer- 
tainty. We can estimate the difference between the two 
ensembles from the inset of Fig. [7)3. At small q around 



Tq the two-point correction is similar in magnitude to ^4, 
but the density fluctuation term becomes negligible for 
wavevectors larger than Qcago- This is consistent with the 
results from facilitated models [5^ . 

The main panel of Fig. [TId shows a prime feature of 
the dynamic susceptibility: its peak height xTil)- 1^ 
corresponds to the maximum in dynamical heterogene- 
ity on a given length scale and thus takes place on times 
of the order of Ta{q). Surprisingly we find xT to have 
appreciable system-size dependence even for a density 
as low as p = 5. At higher densities these effects are 
also pronounced, but their study becomes rapidly com- 
putationally intractable. Transport coefficients analysis 
in Sec. IIIIBI did suggest that a dynamical length scale 
might be as large as the box size for p > 20. But even 
for a system 16 times larger than the typical size con- 
sidered so far and at much lower density, xlil^ t) has not 
yet converged to its bulk value. Considering Xi* does not 
change this observation. Also, not only does xl* keeps in- 
creasing with system size, but it keeps shifting to smaller 
wavevectors. There thus exists a dynamical length scale 
in this system that is much larger than the system size, 
even at densities where caging barely interferes with the 
diffusive regime. Moreover this takes place as the peak 
height, which scales with the dynamical heterogeneity 
volume, remains modest. Such large scale dynamical het- 
erogeneity could result from low-amplitude, long-range 
fiuctuations of the two-point correlation, since their in- 
tegration over a large volume would give them a promi- 
nent contribution. This could then blur the details of 
local dynamical heterogeneity normally associated with 
a dynamical slowdown. Whatever its origin, this effect 
prevents us to quantify completely the wavevector de- 
pendence of X4*((7,i) or the scaling of its peak height, as 
was done in Refs. [Ha, Issj l . A comment remains nonethe- 
less in order. The broad distribution of wavevectors over 
which the peak of xT develops indicates that relaxation 
processes leading to structural relaxation take place over 
a range of length scales. Because no single microscopic 
scale dominates, the mean-field cage opening picture for 
diffusion might be more caricatural than in structural 
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glass formers. Many different microscopic mechanisms 
are probably at play, as the small value of the stretching 
exponent had already suggested in Sec. IIIIBI 

For microscopic q finite-size effects are less important, 
so we will only consider these smaller length scales to 
test theoretical predictions on the other properties of the 
dynamical susceptibility. The full time and wave-vector 
dependence of xlilit) shows a rich structure [53, 55, 58]. 
At short times the motion is ballistic Xiil^t) ~ it 
exhibits a maximum at t*{q) close to the structural re- 
laxation time Ta{q), and at long times goes to unity. Be- 
tween the ballistic regime and the peak the function is 
often fitted to a power-law xlil^ t) ~ since theoreti- 
cal predictions for j{q) differ depending on the dynamical 
relaxation mechanism involved. If short-lived events are 
responsible for the loss of correlations 7 = 1, for indepen- 
dently diffusing defects 7 = 2, while MCT predicts it to 
be the same as the exponent b from the von Schweidler 
form of Eq. [TJ This last scenario is observed in the binary 
L J glass former , but numerical results for kinetically 
constrained models are consistent with the assumption 
of diffusive point-like defects with an anomalous diffu- 
sion exponent [55|. 

Well-separated power-law regimes and the peak of 
xlil^t) can be seen in Fig. [Tj:. The exponent 7, ob- 
tained for wave vectors where the power-law growth lasts 
at least one time decade at p = 20, depends strongly on 
q (Fig. [4}d). To check the MCT prediction we compare 
7 to exponent b extracted from the fit to Eq. [TJ The 
two exponents are significantly different from each other 
for all wavevectors. However, since the von Schweidler 
functional form does not satisfyingly describe the late 
/3 regime even at the highest density considered this is 
not a conclusive assessment. Instead, because of the 
improperly-defined plateau 7 probably corresponds to 
exponent /3 of the stretched-exponential decay, as field- 
theoretic arguments suggest [53] ■ Figure [3Jd presents a 
remarkable agreement between 7 and which support 
this interpretation. A clear separation between the von 
Schweidler and the KWW regimes develops only at den- 
sities higher than what is accessible through simulations, 
so it cannot be excluded that an additional power-law 
regime corresponding to the von Schweidler regime then 
be observed. 



IV. CONCLUDING REMARKS 

We have considered a system of particles formed by 
fixing three orthogonal line segments rigidly at their mid- 
points. Absence of excluded volume implies an absence 
of static correlations, so all the static and thermody- 
namic properties are that of an ideal gas. However, the 
non-crossing condition for the line segments gives rise to 
highly nontrivial dynamics and exhibits "glassy" features 
as the number density is increased. A volume needs to 
open up for a particle to diffuse away from its neighbor 
cage, and this activated dynamics makes the model a 



"strong" glass former. In spite of the inapplicability of 
standard MCT for this system we observe properties that 
are traditionally considered to be success of MCT, such 
as the rescaling of the stretched exponential relaxation in 
Fs{q,t). It remains unclear why such predictions should 
hold here and if some of them break down at densities 
beyond what is computationally reasonable. Note also 
that a model with a similarly trivial static, but fragile 
glass-forming behavior, would also be of great interest to 
test the assumptions that underlie the categorization. 

With increasing density particle displacements acquire 
strong non-Gaussian features on the structural relaxation 
timescale. During this time a small fraction of the par- 
ticles show a much larger MSD than rest. We find these 
"mobile" particles to be associated with local low density 
regions and to cluster. However, the mobile clusters tend 
to be small and highly localized. Yet both the transport 
coefficient decoupling and the system-size dependence of 
the dynamical susceptibility indicate that a sizeable dy- 
namical length scale is present in the system. In light of 
the mobility study and the magnitude of the dynamical 
susceptibility this comes as a surprise, because these are 
usually taken as indirect probes of the dynamical hetero- 
geneity volume. The task to reconcile the large dynami- 
cal length scale with the small size of the mobile regions 
might require identifying a different microscopic metric 
for dynamical heterogeneity. Alternately, this behavior 
shows features that are reminiscent of elastic relaxation 
of a solid after a local volume change. Though this ef- 
fect has not been observed in other glass-forming systems 
so far, it might have been obscured by a stronger local 
dynamical heterogeneity. In any case, a better under- 
standing of this phenomenon would benefit the study of 
all glass-forming systems. 
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APPENDIX: COLLISION FREQUENCY 

We take two needles of length a. In an interval of time 
At, the number of collisions for these two needles is 



Tr,n = 2pvTAt\sm9\a\ 



(A.1) 



where v'^^ is the relative perpendicular velocity and 6 
is the angle between the two Une segments. The factor 
of two appears because two lozenges of size cr^ sin 9 are 
formed. The perpendicular relative velocity averages to 



\ Snrrir 



1/2 



2^' 



(A.2) 



where is the reduced mass and the last equality fol- 
lows from using reduced units. Using (|sin0|) = 7r/4 we 
get 



(A.3) 



Since each cross is made up of three needles, an additional 
factor of 9 has to be included to get the cross collision 
frequency that is used in the text 



-Tec 9 X 1^7171 



9pV7r 



(A.4) 
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